function [y] = m3fSoftmax(xy,theta)
  [n,m] = size(xy);
  [n1,l1] = size(theta);
  if n ~= n1
    error('sizes of xy and theta don''t match');
  end
  y = ones(n,m);
  for i=1:l1
    y = y + (xy >= theta(:,i*ones(1,m))); % modified: Minjie Xu
  end
